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Abstract 

We use A-body simulations to study the evolution of the orbital eccentricities of stars deposited near (< 0.05 
pc) the Milky Way massive black hole (MBH), starting from initial conditions motivated by two competing 
models for their origin: formation in a disk followed by inward migration; and exchange interactions involving 
a binary star. The first model predicts modest eccentricities, lower than those observed in the S-star cluster, 
while the second model predicts higher eccentricities than observed. The A-body simulations include a dense 
cluster of IOMq stellar black holes (SBHs), expected to accumulate near the MBH by mass segregation. Per- 
turbations from the SBHs tend to randomize the stellar orbits, partially erasing the dynamical signatures of 
their origin. The eccentricities of the initially highly eccentric stars evolve, in 20 Myr (the S-star lifespan), to 
a distribution that is consistent at the ~ 95% level with the observed eccentricity distribution. In contrast, the 
eccentricities of the initially more circular orbits fail to evolve to the observed values in 20 Myr, arguing against 
the disk migration scenario. We find that 20%-30% of the S-stars are tidally disrupted by the MBH over their 
lifetimes, and that the S-stars are not likely to be ejected as hypervelocity stars outside the central 0.05 pc by 
close encounters with stellar black holes. 

Subject headings: black hole physics — galaxies: nuclei — stars: kinematics 



1. INTRODUCTION 

In recent years, high resolution observations have revealed 
the existence of many young OB stars in the Galactic cen- 
ter (GC). Accurate measurement of the orbital parameters of 
these stars gives strong evidence for the existence of a mas- 
sive black hole (MBH) which dominates the dynamics in the 
GC (iGhez et al.lll998t HTsenhauer et alJl2005b iGillessen et al l 
2008). Most of the young stars are observed in the central 0.5 
pc around the MBH. The young star population in the inner 
0.05 pc (the "S-stars") consists exclusively of B-stars, in an 
apparently isotropic distribution around the MB H, with rel- 
ativel y high eccentricities (0.3 < e < 0.95; Gilless en et al.l 
2008). The young stars outside this region comprise O-stars 
in one or t wo disks, and present markedly different orbital 
properties (|Levin & Beloborod ov 2003; Ba rtko et al.l 120081: 
iLu et alJl2009F ~ 

Since regular star formation in the region near the MBH is 
inhibited by tidal forces, many suggestions have been made 
regarding the origin of the S-stars. Many of these are proba- 
bly ruled out by observations and/or by the oretical arguments 
(see [Alexan der 2005; Paumard et al. 20061 for a review). The 
various scenarios for the origin of the S-stars predict very dif- 
ferent distributions for their orbits, which in principle could 
be constrained by observations. However, it is not clear to 
what extent relaxation processes can produce changes in the 
distribution of orbital parameters after the stars have been de- 
posited near their current locations. Here we try to resolve 
this question. We use A-body simulations to follow the evo- 
lution of stellar orbits around the GC MBH for 20 Myr, start- 
ing from various initial conditions that were motivated by 
different models for the origin of the S-stars. However, we 
do not discuss here the possibility that an intermediate mass 
black hole was involved in the production and/or evolution 
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of the S-stars, wh ich is discussed in details elsewhere (see 
iMerritt et al.ll2009l and references therein). We find our N- 
body results to be consistent with our analytical predictions, 
and compare them with current observations. We then discuss 
the implications for the validity of the models for the produc- 
tion of the S-stars. In addition we study the possible ejection 
of the S-stars outside the inner 0.05 pc and the contribution 
of such ejected stars to the population o f hypervelocity stars 
(as suggested bv lO'Learv & Loebl d2007| )) and to the isotropic 
population of B-stars observed at distances of up to 0.5 pc 
from the MBH (i.e. at distances similar to those of the young 
O/WR stellar disks, but outside of these disks at high inclina- 
tions). 

In §2 we summarize the different models for the origin of 
the S-stars and the predictions that they make for the stellar 
orbits at the time when the stars are first deposited near the 
MBH. §3 describes the A-body simulations we carried out to 
follow the S-star orbital evolution starting from these initial 
conditions. §4 and §5 present the results of these simulations 
and discusses the implications for the origin of the stellar pop- 
ulations of B-stars in the GC and for the population of hyper- 
velocity stars observed in the Galactic halo. §6 sums up. 

2. MODELS FOR THE S-STARS ORIGIN 

Many solutions have been suggested for the origin of the 
S-stars, but many of these have been e ffectively excluded (see 
lAlexan der 2005it IPaumard et aT1l2006l for a review). Here we 
focus on two basic models which differ substantially in their 
predictions for the initial orbital distribution of the S-stars 
and/or the time passed since their arrival/formation at their 
current location. These are (1) formation of the S-stars in a 
stellar disk close to the MBH, followed by transport through 
a planetary- migration-like scenario to their current positions 
dLevinll2007l) : (2) formation of the S-stars in binaries far from 
the MBH followed by scattering onto the MBH by massive 
perturbers (e.g. gia nt molecular clouds) and tidal disrup- 
tion o f the binaries dPerets et al.l [2007; Perets & Alexanderl 
120071) . leaving a captured star in a tight orbit around the 
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MBH. Binary disruption scenarios similar to (2) have been 
proposed, in which the S-stars formed in a stellar disk (ei- 
ther the currently observed 6 Myr old disk, or an older, 
currently not observed disk) and later changed their orbits 
due t o coherent torques through an instability of eccentric 
disks (Madigan et al. 2008); or through the K ozai mechanism 
result ing from the presence of two disks dLockmann et al.l 
2008). The latter alternative is unlikely since the Kozai mech- 
anism is quenched in the presen ce of a massive enough cusp 
of sta rs such as exists in the GC (Chang 2008; Madig an et al.l 
2008)). In any case, all the binary-disruption scenarios imply 
very similar initial distributions for the captured S-stars, and 
may differ only in their relevant timescales. 

In the following, we briefly discuss the initial distribution 
of the eccentricities and inclinations of the S-stars expected 
from the different scenarios for their production. The models 
are summarized in Table Q] 

2.1. Binary disruption by a massive black hole 

A close pass of a binary star near a massive black hole re- 
sults in an exchange interaction, in which one star is ejected 
at high velocity, while its companion is captured by the MBH 
and left on a bound orbit. Such interaction occurs because 
of the tidal forces exerted by the MBH on the binary compo- 
nents. Typically, a binary (with mass, Mu n , and semi-major 
axis, abin), is disrupted when it crosses the tidal radius of 
the MBH (of mass M»), given by rt — abiniMt/Mbin) 1 ^ ■ 
One of the binary c omponents is captured by the MBH 
(iGould & Ouillenl l2003 ) on a wide and eccentric orbit while 
the companion is ejected with high velocity (Hills 1988|). 

The capture probability and the semi-major axis distribu- 
tion of the captured stars were estimated by means of nu- 
merical simulations, showing that most binaries approa ching 
the MBH within the tidal radius rtjabin) are disrupted (Hills 



fT99llfT99llGualandris et al.l2005tlBromiev et al.l2006l) . The 
harmonic mean semi-major axis for 3 -body exc hanges with 
equal mass binaries was found to be (iHillsll 19911) 
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(i) 

where aun is the semi-major axis of the infalling binary and 
a cap that of the captured star (the MBH-star "binary"). Most 
values of a cap fall within a factor 2 of the mean. This relation 
maps the semi-major axis distribution of the infalling binaries 
to that of the captured stars: the harder the binaries, the more 
tightly bound the captured stars. The periapse of the cap tured 
star is at r t , and therefo re its eccentricity is very high riffiUsl 
[19911 iMiller eTaT]l200lh . 

e=l - r t /a cap ~ 1 - 1.8(Af bin /Af.) 1 / 3 ~ 0.94 - 0.99 (2) 

for values typical of B-type main sequence binaries and the 
MBH in the GC (M bm = 6-30 M ; M. = 3.6 x 10 6 M ; 
(leap = 0.5 — 2 x (a cap )). Therefore, in order to study the 
evolution of S-stars from the binary disruption scenarios we 
assume that the initial eccentricities of S-stars are in the range 
0.94 — 0.99 (where most are close to the mean value of 0.98). 

In principle the binary disruption scenario has specific pre- 
dictions for the semi-major axis distribution of the captured 
stars, which could also be used for constraining the model. 
However such distribution is highly sensitive to differences in 
the (unknown) binary distribution in the GC region. The pre- 
diction of high eccentricities for the captured S-stars, instead, 



is robust and has only a weak dependence on the mass of the 
binary. 

The inclinatio ns of the captured S-stars in the massive per- 
turbers scenario (Per ets et al.l2007l) are likely to be distributed 
isotropically since the stars originate in an isotropic cusp. Al- 
though in the disk instability scenario (Madi gan et aT1l2008l) 
the progenitors of the captured stars form in the stellar disk, 
their original inclinations could be excited to higher inclina- 
tions than the typical inclinations observed in the stellar disk 
(Y. Levin, private communication), and may resemble a more 
isotropic distribution. 

2.2. Planetary like migration from the young stellar disk(s) 

iLevinl (120071) suggested that the S-stars could have formed 
in the currently obser ved stellar disk in the GC (Bart ko et al.l 
2008] iLuet alJl2009h . and then migrated inward in a way 
similar to planetary migration. The migration time scale ex- 
pected from such a scenario could be as short as 10 5 yr (for 
type I migration), w hich could be compara ble to (although 
possibly larger than; Nayaksh in et al.l (l2007|)) the lifetime of 
the ga s eous disk. Recen t analy tic work (lOgilvie & Lubowl 
d2003l) : iGoldreich & Saril d2003l) ; and references therein) has 
shown that eccentricity is likely to be damped during migra- 
tion, unless eccentricity excitation occurs, which requires the 
opening of a clean gap in the dis k. In the latt er case the migra- 
tion timescale might be larger (lLevinll2007l), possibly incon- 
sisten t with the lifetime of the gaseous disk ( Nayakshin et al.1 
120071) . It is therefore more likely that the eccentricities of 
the stars are damped during the migration. Even eccentricity 
excitation, if such took place, is unlikely to excite very high 
eccentricities. The mean eccentricity of the observed stars in 
the stellar disk is 0.34 ± 0.06. Therefore, in order to study 
the evolution of the S-stars following their formation in a stel- 
lar disk, we assume them to have low eccentricities, or, be- 
ing conservative, moderately high eccentricities (e max = 0.5; 
where we use a thermal distribution of eccentriciti es, cut off at 
emax)- These simulations include the less likely dLevin 2007) 
possibility that the stellar disk extended inward to the current 
region of the S-stars, in which case the S-stars were formed 
in-situ and did not migrate. 

3. THE N-BODY SIMULATIONS 

To test these competing models, we carried out A-body 
simulations of the inner Milky Way bulge using mod- 
els containing a realistic number of stars. All integra- 
tions were carried out on the 32-node GRAPE cluster 
gravitySimulator at the Rochester Institute of Tech- 
nology which adopts a parallel setup of GRAPE accelerator 
boards to efficiently compute gravitat ional forces. The d irect- 
summation code </>GRAPE was used (lHarfst et al.l l2007). The 
simulations used a softening radius of ~ 4 Rq comparable to 
the radius of the S-stars r*, so as to be able to follow even the 
closest encounters between stars. 

Our initial conditions were based on a collisionally evolved 
model of a cusp of stars and stellar rem nants around the 
GC MB H dHopman & Alexandeil l2006bl hereafter HA06; 
see also Freita get al.ll2006l) . HA06 evolved the multi-mass 
isotropic Fokker-Planck equation representing the stellar dis- 
tribution in the region extending to ~ 1 pc from the MBH; 
in their models the contribution to the gravitational potential 
from the distributed mass is ignored. HA06 fixed the relative 
numbers of objects in each of four mass bins by assuming a 
mass function consistent with continuous star formation. The 
10M Q stellar-mass black holes (SBHs) were found to follow 
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Table 1 

Models for the S-stars origin 



# 


Origin 


Initial 


Time 


Model 


Survival 


Refs. c 






Eccentricity 


(Myr) Probability" Fraction 6 




1 


Capture following binary disruption due 


High 


6 


0.26 


0.8 


1 




to disk instability in the currently observed disk (0.94 < e < 0.99) 










2 


Capture following binary disruption due 


High 


20 


0.93 


0.7 


2 




to massive perturbers or to disk instability 


(0.94 < e < 0.99) 












in an old non-observed disk 












3 


Disk formation + planetary like 


Low 


6 


8 x 10" 3 


0.9 


3 




migration (currently observed disk) 


(e < 0.5) 










4 


Disk formation + planetary like 


Low 


20 


0.06 


0.8 


3 




migration (possible old disk) 


(e < 0.5) 











"Probability for the samples of the observed and simulated S-stars 

to be randomly chosen from the same distribution (see text). 

^Fraction of S-stars not disrupted by the MBH during the simulation (see text). 

c (1) Madigan et al. (2008) (2) Perets et al. (2007) (3)Levin (2007) 



a steep, n(r) ~ r~ 2 density profile near the MBH while the 
lower-mass populations (main-sequence stars, white dwarfs, 
neutron stars) had n ~ r~ a , 1.4 < a < 1.5. The SBHs were 
found to dominate the mass density inside ~ 0.01 pc. 

Based on this model, we constructed an iV-body realization 
containing a total of 1200 objects within 0.3 pc of the MBH: 
200 "stars," with masses of 3Mq, and 1000 "black holes" 
with masses 10 M , around a MBH of 3 x 10 6 M . We set 
a = 2 for the SBHs and a = 1.5 for the lower-mass stars, and 
each density component was tapered smoothly to zero beyond 
0.1 pc when computing the corresponding f(E). Since the S- 
stars may have masses as high as ~ 10 Mq, the higher mass 
stars in our simulations could also be treated as S-stars. We 
did not see any major differences in the evolution of the more 
massive and the less massive stars, and we discuss the evolu- 
tion of both together. 

The number of SBHs contained within a radius r in our N- 
body models was 

iV(< r) » 600 ( 1 r^—) (3) 

Vo.i pcy 

implying a distributed mass within 0.1 pc of ~ 10 4 M Q . This 
is somewhat (~ 2 — 3x) lower than the mass in SBHs in 
the HA06 or simila r (Morris 1 993t iMiralda-Escude & Gould! 
2000; Freitag_etaL|2006) models at the same radius, and a 
factor ~ 5 lower than the total mass (mostly in main-sequence 
stars) in the Fokker-Planck models. In this sense, the rates of 
evolution that we infer below can be considered to be conser- 
vative. 

On the other hand, we note that the late-type (old) stars that 
dominate the number counts in this region have a much flatter 
density profile than predicted by the HA06 models, possibly 
even exhibiting a central "hole" dFiger et aLll2003t IZhu et all 
2008). Only the B-type stars in the nuclear sta r cluster show 
a steeply-rising number de nsity, a = 1.1 ±0.3 dSchodel et alj 
l2007tlGillessen et alj|2008l) but they presumably constitute a 
negligible fraction of the total mass in this region, and in any 
case are far too young to have reached a collisional steady 
state around the MBH. While the origin of this discrepancy 
between models and observations is currently unresolved, it 
may imply that the other contributors to the distributed mass 
around the MBH, including the SBHs, also have a lower den- 
sity than in the Fokker-Planck models. For instance, relax- 
ation times at the GC may be too long for collisionally re- 
laxed steady states to have been established in the last 10 Gyr 
dMerritt & Szellll2006l) . 



Modeling of the stellar proper motion data (Trip pe~et alj 
l2008t ISchoedel et al.ll2009h implies a distributed mass within 
1 pc of 0.5 — 1.5 x 10 6 Af Q , but these data are consistent with 
both rising and falling mass densities within this region and 
the distrib uted mass in the inn er 0.1 pc is essentially uncon- 
strained (ISchoedel et alJl2009t) . 

Because of these uncertainties, we discuss below how our 
results would vary if different numbers of SBHs were as- 
sumed. 

As discussed in the previous section, we studied two basic 
sets of initial conditions for the S-stars. In the first model 
we assumed that the S-stars w ere captured by the MB H as in 
the binary disruption scenario dGould & Ouillenl2003h . which 
leaves the captured stars in highly eccentric orbits (> 0.94 — 
0.99; cf. section |2~TT >. Under these assumptions, the stars 
evolved for 6 Myr (if formed in the stellar disk) or longer (if 
the S-stars formed outside the central pc, i.e. not in the young 
stellar disk). We evolved the models for up to 20 Myr, which 
is comparable to the lifetime of the observed S-stars (although 
some may have longer lifespans). In the second scenario we 
assumed that the S-stars formed in a gaseous disk and then 
migrated inwards (or formed in situ in a disk extending close 
to the MBH). For this case we assumed the S-stars to have low 
eccentricities (< 0.5), as typical of disk formation models, 
and to evolve for 6 Myr (the lifetime of the observed stellar 
disk). In order to check both scenarios we selected the stars 
with initially high eccentricity orbits (0.94 < e < 0.99) and 
low eccentricity orbits (e < 0.5) and followed their evolution 
for a time appropriate to their presumed origin. 

In addition to these evolutionary scenarios we also studied 
the possibility of ejection of S-stars as hypervelocity stars, 
following a close encoun ter with a SBH in the vicinity of 
the MBH, as suggested by lO'Leary & Loebl (120071) . Our high 
resolution simulations can accurately follow close encounters 
between stars, and therefore track any resulting high velocity 
ejections of stars. Motivated by recent observations of B-type 
main sequence stars outside the inner 0.05 pc of the Galaxy, 
which show evidence of random orientations similar to those 
of the S-stars, we also looked for stars scattered to larger dis- 
tances by gravitational encounters. In other words, we con- 
sidered whether captured S-stars could dynamically evolve to 
become the extended B-type stars population observed out- 
side the central 0.05 pc. 

Although some binaries cou ld exist in the close regions near 
the MBH (lPeretsll2008l l2009h . these are likely to be rapidly 
disrupted and not play an important role in the dynamical 
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evolution of the S-stars (Perets 2009). We therefore do not 
include any binaries in our initial conditions for the S-stars 
and SBHs evolution. 

4. RESULTS 

4.1. Simulations vs. theory: resonant relaxation 
We applied the correlation curve method lEilon et all 



(2008) to our simulations to identify the relaxation pro 
cess responsible for the dynamical evolution of the stars 
(Fig. [TJi. The method is able to detect and measure 
relaxation in nearly Keplerian A-body systems. In the 
isotropic system considered here, the angular momentum 
of the stars, J, evolves bo th due to the slow stoc hastic 
two-body relaxation (e.g. iBinney & Tremainel 1 19871) and 
to the rapid resonant relaxa t ion (RR) dRauch & Tremainel 



1996; Rauch & Ingalls 1998; Hopman & Alexander 2006a; 
Giirkan & Hopm an 2007: lEilon et al.1 12008). Two-body re- 
laxation changes J in a random walk fashion, |AJ| / J c — 
\Jt/tnr over the long two-body relaxation time-scale 
t~n b ^Q 2 / (N log Q), where J c is the maximal (circular or- 
bit) angular momentum for a given energy, Q — Af./M*, 
N is the number of enclosed stars on the distance scale of 
interest, and time r = t/P is measured in terms of the or- 
bital period on that scale. Resonant relaxation occurs when 
the symmetries of the potential act to constrain the stellar 
orbits (e.g. closed ellipses in a Kepler potential, or planar 
rosettes in a spherical one). As long as the symmetry is ap- 
proximately maintained on the coherence timescale t w , the 
stars experience coherent torques, and AJ/J C ~ (VN /Q)t. 
In a nearly Keplerian potential, as is the case in the inner 
parsec of the GC, RR can change both the magnitude of J 
("scalar RR") and its direction ("vector RR"). In the New- 
tonian context, the coherence of scalar RR is limited by the 
precession of the apoapse due to the enclosed stellar mass, 
on a time-scale tm ~ Q/N. On timescales r > tm, the 
coherent change A J(tm) becomes the mean free path in J- 
space for a rapid random walk, \AJ\ / J c = yr~prsRR, where 
t 8 rr ~ Q <§C tjv r- The coherence of vector RR is self-limited 
by the change in the orbital orientation due to RR, and is even 
faster, | AJJ / J c = ^Jt/t vRR , where t vRR ~ Q/-/N. 

Figure (j} shows the rms change in the scalar and vector 
angular momentum of the stars in the simulation, as a func- 
tion of the time lag r (correlation curves), up to the maxi- 
mal time lag for which the simulation can still be analyzed 
with high statistical confidence, r max ~ 10 4 . In our simula- 
tion 3 , the scalar RR coherence time expected from theory is 
tm ~ 283 and the scalar RR timescale is t sRR ~ 3.1 x 10 5 . 
The behavior of the scalar correlation curve clearly indicates 
that RR dominates relaxation in our GC model. The turn from 
the coherent phase to the random- walk phase of RR is seen at 
r ~ 300, close to the predicted value of tm- The random- 
walk growth continues up to full randomiz ation and satura- 
tion at ~ 70, as expected (Eil on et aDl20 08). The scalar RR 
timescale can be estimated directly from the simulated data 
by trr = r/(|AJ| /J c ) 2 ; substituting \AJ\/J C ~ 0.15 at 
r 5000 yields t s rr~ 2.2 x 10 5 , in good agreement 

with the predicted value. 

Furthermore, since J/J c = ^/T~~i?, the RR-driven eccen- 
tricity evolution relates the final eccentricity, e/, after time 

3 For a spectrum of masses, tm ~ Q/N (M») and t 3 rr ~ M»/M e g, 
where (M») =8.8 M© and M cff = {M%) / (M»). M eff ~9.7M© for our 
simulation. 
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FIG. 1 . — The correlation curves | AJ\ / J c and | AJ| / J c as a function of 
the time lag r for all the stars in the simulation. The simulation data (points) 
are compared to t he best fit theoretical curves (lines). Four regimes are seen 
lEilon et al. 2008, see text and also). At very short timescales, non-resonant 
two-body relaxation (<x y/r) dominates over RR (<x t); At r < tm the curves 
display a coherent, linear rise; At tm < r < ta, the curves display a random- 
walk growth and at r > tj, they begin to saturate. 



lag t to the initial eccentricity, e;. We therefore expect the 
eccentricities at some given time lag to have a typical spread, 



e ( r ) ^5 e /( r ) ^ e+ ( r )' where 



T 

trr 



1/2 



(4) 



The magnitude of the predicted change in eccentricity agrees 
well with that observed in the simulations. For example, an S- 
star initially captured by a tidal exchange event on a P = 500 
yr (a ~ 0.04 pc), e, = 0.97 orbit can evolve by RR to an 
e/ = 0.80 orbit in 20 Myr. The short vector RR timescale 
t v rr ~ 10 4 (t vRR — 5x 10 6 yr at 0.04 pc) implies full ran- 
domization of the orbital planes after 6 Myr, throughout the 
S-cluster volume, as is observed in the simulation. We there- 
fore conclude that RR is the dominant mechanism responsible 
for the dynamical evolution of the S-stars and other stars close 
to the MBH in the GC. 

4.2. S-star eccentricities and inclinations 

In Figure [2] we show the final cumulative eccentricity dis- 
tribution of the S-stars for the different origin models (Ta- 
ble [TJ. These are co mpared to the the or bits of the observed 
S-stars (taken from Gillessen et al. 2008). The probabilities 
for the samples of the observed and simulated S-stars to be 
randomly chosen from the same distribution (calculated us- 
ing a two-sample \ 2 test ) &g given in Table [TJ We find that 
the binary disruption model taking place at least 20 Myr ago 
is much favored over all other models tested here. We find 
93% chance for the observed S-stars and the simulated stars 
in such model to originate from the same distribution (with 
the shorter timescale of 6 Myr consistent at the 26% level) . 
In contrast, the disk migration scenarios seem to be excluded 
(for the given assumptions), since they have major difficulties 
in explaining the large fraction of eccentric orbits observed 
for the S-stars in the GC. 

We find that the inclination distribution of the cap- 
tured stars, initialized with the same inclination, is rapidly 
isotropized, to resemble a random distribution of inclinations 
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— Observed S-stars (Gillessen et al. 2008) 
Simulated S-stars (e init <0.5) after 6 Myrs 

Simulated S-stars (e. nit <0.5) after 20 Myrs 

- - -Simulated S-stars (0.94<e<0.99) after 6 Myrs 
Simulated S-stars (0.94<e. <0.99) after 20 Myrs 



■ Thermal eccentricity distribution 
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FIG. 2. — Cumulative distribution of observed and simulated S-stars eccen- 
tricities, for various models (see legend). 



(consistent with random at the ~ 65% and ~ 20% level, af- 
ter evolution of 20 and 6 Myr, respectively). This is expected 
from the RR process, as discussed in the previous section (see 
also fig. [T). We conclude that the observed isotropic distribu- 
tion of the S-stars angular momentum direction is consistent 
with all the S-stars production models studied here, and can 
not be used to discriminate between them, although it con- 
strains the lifetime of the S-stars system to be at least ~ 4 
Myr at the 95% level, assuming all S-stars were initially put 
on the same plane. 

4.3. Survival of the S-stars: tidal disruption, ejection and 
hypervelocity stars 

As already discussed above, the S-stars can change their 
orbits due to their dynamical evolution. A star could there- 
fore be scattered very close to the MBH and be disrupted by 
it, if its pericenter distance from the MBH becomes smaller 
than the tidal radius of the star r t = r^M./m*) 1 / 3 . Many 
of the S-stars could therefore not survive for long close to 
the MBH. We followed the orbits of stars in our simulations 
and calculated the fraction of stars that have been disrupted. 
For the tidal disruption calculations all stars were assumed 
to have the typical main sequence radius according to their 
mass. We consider a star as being disrupted if its pericenter 
became smaller than twice the tidal radius during the simula- 
tion (i.e., when it is strongly affected by the MBH tidal forces 
or even totally disrupted in one pericenter passage). We find 
that most of the S-stars survived to current times in all the 
models (see survival fractions inHJ. The S-stars population in 
the GC therefore gives a good representation of all the S-stars 
formed/captured in this region. The production rate of the S- 
stars required to explain current observations is therefore only 
slightly higher (1.2— 1.3 times higher) than that deduced from 
current number of S-stars observed. 

In principle the S-stars could be ejected by strong 
encounters to orbits with larger semi-major axes, 
putting them outside the 0.05 pc region near the MBH 
( Miralda-Escud 6 & Gouldl l2000|), or even ejecting th em as 
unbound hypervelocity stars (lO'Leary & Loebl 12007). The 
softening radius used in our simulation was r so f t = 4R & , 
comparable to the radius of observed S-stars, allowing us to 



follow even very close encounters. Nevertheless, we find that 
only ~ 10% of the stars were ejected outside of the central 
0.05 pc, and even those had maximal semi-major axis in 
the end of the simulation not extending beyond 0. 1 pc. We 
therefore conclude that such ejected S-stars can not explain 
recent observations of many B-type stars outside the central 
0.05 pc (H Bartko, private communication). Moreover, none 
of the 3 Mq stars in our simulations have been ejected as a 
hypervelocity star, suggesting that ejection of hypervelocity 
stars through encou nters with S BHs is not an efficient 
mechanism (see also Perets 2009 for the related constraints 
on this mechanism). We note that since our simulations 
can be rescaled, we can probe much higher stellar densities 
and smaller softening radius (see next section). When such 
rescaling is used, in which all the 1200 stars are distributed 
between 3 x 10~ 4 pc to 0.1 pc (rescaling the radii by half), we 
find 5 stars ejected beyond 0.1 pc (to have final semi-major 
axis of 0.1, 0.16, 0.16, 0.21 and 0.37 pc), but none ejected 
as hypervelocity stars or even as slow unbound stars. 

5. DEPENDENCE OF THE RESULTS ON THE ASSUMED DENSITY 
OF SBHS 

As discussed above, the density of the SBHs that are re- 
sponsible for the evolution in our A-body models is not well 
determined. Scaling the A-body results to different assumed 
values of A is complicated by the fact that scalar resonant re- 
laxation has two regimes, coherent and random walk. In our 
models, the transition occurs at 

where the Aw H scaling is for scattering by SBHs of a given 
mass. 

In the coherent regime, At < tu, orbital angular momenta 
grow as 



AJ 4 (N SBH \ 1/2 At 

J c ~ \ 10 3 J P 
In the diffusive regime, At > tM, 



AJ 



Jc 



2 x 10 _3 W — 



(6) 



(7) 



T S R.R 

independent of Asbh- 

We are interested in the orbital evolution of the S-stars 
over timescales of At ~ O(10 7 ) yr. In our simulations and 
those with larger N, At ijw. In this large- A regime, 
changes in eccentricity are dominated by the diffusive rela- 
tion and are therefore expected to be nearly independent of A 
for time scales of interest, at least up to A-values of ~ 10 5 
where the distributed mass begins to approach the mass of the 
MBH and resonant relaxation is no longer effective. Only for 
AgBH ^ 10 does tjvf approach 10 7 yr and our results start de- 
pending significantly on Asbh- However, such a small num- 
ber of SBHs in the volume of interest is highly unlikely, and 
therefore our results are robust to the details of the SBH cusp 
model. 

The A-body simulations can also trivially be rescaled by 
r^Ar, t-> A 3/2 t (8) 
at fixed mass. This corresponds to placing the same number 
of SBHs into a smaller (larger) region and integrating for a 
shorter (longer) time. For instance, if we rescale our simula- 
tions to three times smaller distances (in which case the stars 
are distributed between 3 x 10~ 4 and 0.05 pc), the integration 
time becomes ~ 5 Myr. 
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6. SUMMARY 

The approximately thermal eccentricity distribution of the 
S-stars near the massive black hole (MBH) in the galactic cen- 
ter (GC), N(< e) oc e 2 , is not naturally predicted by either 
of the two leading models for their production: migration of 
stars formed in a gaseous disk; or capture of stars following 
binary disruption by the MBH. The former model predicts ec- 
centricities that are too low, the latter too high. In this paper, 
we followed the dynamical evolution of orbits of various ec- 
centricities near the GC MBH, including for the first time the 
cluster of stellar-mass black holes (SBHs) that is expected to 
form around the MBH via mass segregation. We found that 
perturbations from the SBHs can reduce the eccentricities of 
initially very eccentric orbits (0.94 < < 0.99) into a dis- 
tribution that is consistent with the observed one, in a time of 
approximately 20 Myr, comparable to S-star lifespans; some 
of the stars change their eccentricities by more than 0.5 to val- 
ues as low as e fi na i — 0.2 — 0.4. We confirmed these TV -body 
results via a theoretical analysis of the relaxation process, and 
used that analysis to argue that our results are not strongly de- 
pendent on the (unknown) normalization of the SBH density 
near the MBH. The same mechanism is unable to convert ini- 
tially low-eccentricity orbits into very eccentric ones on the 



same time scale, arguing against the validity of the disk mi- 
gration model for the origin of the S-stars. We also found that 
most S-stars are not disrupted by the MBH during their life- 
time, and very few are ejected outside the central 0.05 pc near 
the MBH, and none having a semi-major axis beyond 0.1 pc. 
We did not find any hypervelocity star ejected in our simula- 
tion. 

Evolution toward a thermal eccentricity distribution is a nat- 
ural consequence of random gravitational encounters with a 
population of massive perturbers. In this paper, we considered 
the effect of a background of SBHs, which are expected on the 
basis of very general arguments to contribute a total mass of 
~ 10 4 Mq in the inner 0.1 pc around the MBH. We showed 
that they were effective at moderating the eccentricities of ini- 
tially highly eccentric orbits. These results strengthen models 
in which the S-stars form from disrupted binaries, and disfa- 
vor models in which the S-stars are formed with low eccen- 
tricities. 
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